function results = equilibrium_price_change(prefix_pf,prefix_skill,prefix,drop_stop,knitro_path)
    addpath('../')
    addpath('../Densities')
    addpath('../UtilityFunctions')
    addpath('../ProductionFunctions')
    addpath('../Economies/')
    addpath(knitro_path)

    dir_in = '../../../data/generated/matlab/Calibration/'
    dir_out = '../../../data/generated/matlab/simulation/'

    w0_range = 1e-4;
    wf_range = 1e4;

    % load calibrated parameters for skill distribution
    par_cal = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_sol.csv']);
    par_mom = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_mom.csv']);
    inequ_par = tools.csv2struct([dir_in,'inequ_par.csv']);
    constants = tools.csv2struct('../../../data/raw/constants.csv');

    DensSkill = abilityDistFixed(par_cal);

    DensPart = ParticipationNormal(par_cal);

    % utility function object
    par_util.epsilon = constants.epsilon;
    Util = UtilityQuasiLinear(par_util);

    % combine in SkillAndParticipation object
    par_dist.dens_skill = DensSkill;
    par_dist.dens_part = DensPart;
    par_dist.util = Util;

    par_dist.num_nodes_phi = 50;
    par_dist.scale_inc = par_cal.scale_inc; 

    par_dist.hsv_tax = true;
    Dist = SkillAndParticipation(par_dist);

    % Production function
    par_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_sol.csv'])
    mom_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_mom.csv'])

    par_pf.parse_inputs = true;
    par_pf.KrusellRobotsEach_use_properties = true;
    ProdFun = KrusellRobotsEach(par_pf);

    % Economy object
    par_econ.pf = ProdFun;
    par_econ.util = Util;

    par_econ.w_lb = 1e-6;
    par_econ.w_ub = 1e4;

    % robot prices and taxes
    inp.tau_B = par_pf.tau_B;
    inp.tau_E = par_pf.tau_E;
    inp.tau_S = par_pf.tau_S;

    par_econ.q_B = mom_pf.q_B;
    par_econ.q_E = mom_pf.q_E;
    par_econ.q_S = mom_pf.q_S;

    q_B_0 = par_econ.q_B;
    q_B_1 = mom_pf.Y_B_M_1/(1+inp.tau_B);

    par_econ.dist = Dist;

    par_econ.solver = 'knitro_eqsolve'

    par_econ.phi_lb = constants.phi_lb;
    par_econ.r = inequ_par.inequ_par;
    par_econ.par_dist = par_cal;

    par_econ.varnames = {'L_M','L_R','L_C', ...
                        'K_B_M','K_B_R','K_B_C', 'K_E','K_S',...
                        'tau_B','tau_E','tau_S','transfer','tau_hsv','lambda_hsv','compensation'};

    econ = KrusellWithRobotsEachEconomy(par_econ);
    econ.set_num_nodes(100,50,20,50);

    inp.L_M = par_pf.L_M_0;
    inp.L_R = par_pf.L_R_0;
    inp.L_C = par_pf.L_C_0;
    inp.K_B_M = par_pf.K_B_M_0;
    inp.K_B_R = par_pf.K_B_R_0;
    inp.K_B_C = par_pf.K_B_C_0;
    inp.K_E = par_pf.K_E_0;
    inp.K_S = par_pf.K_S_0;
    inp.transfer = par_cal.transfer;
    inp.tau_hsv = par_cal.tau_hsv;
    inp.lambda_hsv = par_cal.lambda_hsv;
    inp.compensation = 0;

    lb.L_M = 1e-16;
    lb.L_R = 1e-16;
    lb.L_C = 1e-16;
    lb.K_B_M = 1e-16;
    lb.K_B_R = 1e-16;
    lb.K_B_C = 1e-16;
    lb.K_E = 1e-16;
    lb.K_S = 1e-16;
    lb.tau_E = inp.tau_E;
    lb.tau_S = inp.tau_S;
    lb.transfer = inp.transfer;
    lb.tau_hsv = inp.tau_hsv;
    lb.lambda_hsv = inp.lambda_hsv;
    lb.compensation = 0;


    ub.L_M = 1e16;
    ub.L_R = 1e16;
    ub.L_C = 1e16;
    ub.K_B_M = 1e16;
    ub.K_B_R = 1e16;
    ub.K_B_C = 1e16;
    ub.K_E = 1e16;
    ub.K_S = 1e16;
    ub.tau_E = inp.tau_E;
    ub.tau_S = inp.tau_S;
    ub.transfer = inp.transfer;
    ub.tau_hsv = inp.tau_hsv;
    ub.lambda_hsv = inp.lambda_hsv;
    ub.compensation = 0;

    drop_range = sort([0:0.001:drop_stop]);

    %% vesion with flexible capital
    % drop_range = 0;
    inp.tau_B = par_pf.tau_B;
    lb.tau_B = inp.tau_B;
    ub.tau_B = inp.tau_B;

    econ.fix_other_capital = false;
    simulate(econ, inp, lb, ub, q_B_0, drop_range,w0_range,wf_range,prefix, 'simulation',dir_out,par_cal)

    %% vesion with flexible capital and high robot tax

    inp.tau_B = 10 * par_pf.tau_B; % 10 times higher tax on robots
    lb.tau_B = inp.tau_B;
    ub.tau_B = inp.tau_B;

    econ.fix_other_capital = false;
    simulate(econ, inp, lb, ub,  q_B_0, drop_range,w0_range, wf_range, prefix, 'simulation_high_tax',dir_out,par_cal);

    %% vesion with fixed capital (as sanity check for the wage-decile elasticities)
    inp.tau_B = par_pf.tau_B;
    lb.tau_B = inp.tau_B;
    ub.tau_B = inp.tau_B;

    drop_range = sort(linspace(0,1-q_B_1/q_B_0,10));
    econ.fix_other_capital = true;
    [results, results_percentiles] = simulate(econ, inp, lb, ub, q_B_0,drop_range,w0_range,wf_range,prefix,'simulation_fix_capital',dir_out,par_cal);
    percentiles_0 = results_percentiles(results_percentiles.q_B==q_B_0,:).value;
    percentiles_1 = results_percentiles(abs(results_percentiles.q_B-q_B_1)<1e-6,:).value;

    percentiles_change = percentiles_1 - ...
        percentiles_0;

    K_B_0 = results(results.q_B==q_B_0,:).K_B;
    K_B_1 = results(abs(results.q_B-q_B_1)<1e-6,:).K_B;

    K_B_rel_change = (K_B_1- K_B_0)/K_B_0;

    dir.generated = '../../../../data/generated/r/';
    dec_eps = tools.csv2struct([dir.generated,'dec_eps_AR2018.csv']);
    dec_eps = dec_eps.eps;

    targets.eps_w_K.val = dec_eps;

    bar([targets.eps_w_K.val,percentiles_change./ ...
         percentiles_0./K_B_rel_change])
    legend('data','model')
end

function [results, results_percentiles] = simulate(econ,inp,lb,ub, q_B_0,droprange, w0_range,wf_range,prefix, suffix, dir_out, par)
    sol = inp;
    stats = struct();
    percentiles = struct();
    results = table();
    results_percentiles = table();

    for i = 1:length(droprange)
        % change price of robots and recompute
        drop = droprange(i);
        econ.q_B = q_B_0 * (1-drop);
        
        fprintf('\n=================================================\n')
        fprintf('%s, price is %f\n', suffix, econ.q_B)
        fprintf('=================================================\n')
        
        if econ.fix_other_capital
            lb.K_E = inp.K_E;
            lb.K_S = inp.K_S;
            
            ub.K_E = inp.K_E;
            ub.K_S = inp.K_S;
        end
        
        tic
        [sol,fval,exitflag,output] = econ.compute_equilibrium(sol,lb, ...
                                                          ub);
        toc
        
        [stats, percentiles] = compute_stats(econ, sol, par);
        stats.pricedrop = drop;
        stats.exitflag = exitflag;
        stats.fval = fval;
        results = [results; struct2table(stats)];
        results_percentiles = [results_percentiles; struct2table(percentiles)];
    end

    results_struct = table2struct(results);
    results_percentiles_struct = table2struct(results_percentiles);
    tools.struct2csv(results_struct,[dir_out,prefix,'_',suffix,'.csv'])
    tools.struct2csv(results_percentiles_struct,[dir_out,prefix,'_',suffix,'_percentiles.csv'])
end

function [stats, percentiles_strct] = compute_stats(econ, sol, par)
    stats = sol;
    stats.K_B = sol.K_B_M + sol.K_B_R + sol.K_B_C;
    stats.q_B = econ.q_B;
    stats.q_E = econ.q_E;
    stats.q_S = econ.q_S;

    [stats.Y_M,stats.Y_R,stats.Y_C,stats.Y_B_M,stats.Y_B_R,stats.Y_B_C,stats.Y_E,stats.Y_S]=econ.factor_prices(sol,par);
    stats.Y = econ.pf.Y(sol,par);

    w_lb = 1e-6;
    w_ub = 1e4;

    stats.part = econ.dist.mass_participants(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,par);
    stats.emp_M = econ.dist.emp_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_R = econ.dist.emp_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_C = econ.dist.emp_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    stats.avg_wage_M = econ.dist.average_wage_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_R = econ.dist.average_wage_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_C = econ.dist.average_wage_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    % welfare
    derived_vars = econ.compute_vars_from_input(sol);
    stats.welfare_part_M = econ.welfare_participants_M(derived_vars);
    stats.welfare_part_R = econ.welfare_participants_R(derived_vars);
    stats.welfare_part_C = econ.welfare_participants_C(derived_vars);

    stats.welfare_non_part = econ.welfare_non_participants(derived_vars);
    stats.welfare_part = econ.welfare_participants(derived_vars);
    stats.welfare = econ.welfare(derived_vars);

    % welfare weights
    stats.tot_welfare_weight_part_M = econ.tot_welfare_weight_participants_M(derived_vars);
    stats.tot_welfare_weight_part_R = econ.tot_welfare_weight_participants_R(derived_vars);
    stats.tot_welfare_weight_part_C = econ.tot_welfare_weight_participants_C(derived_vars);
    stats.tot_welfare_weight_nonpart = econ.tot_welfare_weight_non_participants(derived_vars);
    stats.tot_welfare_weight_part = ...
        stats.tot_welfare_weight_part_M + ....
        stats.tot_welfare_weight_part_R + ...
        stats.tot_welfare_weight_part_C;
    
    stats.tot_welfare_weight = ...
        stats.tot_welfare_weight_part + ...
        stats.tot_welfare_weight_nonpart;
    
    
    stats.avg_welfare_weight_part_M = stats.tot_welfare_weight_part_M/...
        (stats.tot_welfare_weight * stats.emp_M * stats.part);
 
    stats.avg_welfare_weight_part_R = stats.tot_welfare_weight_part_R/...
        (stats.tot_welfare_weight * stats.emp_R * stats.part);

    stats.avg_welfare_weight_part_C = stats.tot_welfare_weight_part_C/...
        (stats.tot_welfare_weight * stats.emp_C * stats.part);   
    
    stats.avg_welfare_weight_nonpart = stats.tot_welfare_weight_nonpart/...
        (stats.tot_welfare_weight * (1-stats.part));
    
    % normalized welfare weights for participants
    stats.avg_welfare_weight_part_norm_M = stats.tot_welfare_weight_part_M/...
    (stats.tot_welfare_weight_part * stats.emp_M);
   
    stats.avg_welfare_weight_part_norm_R = stats.tot_welfare_weight_part_R/...
    (stats.tot_welfare_weight_part * stats.emp_R);
    
    stats.avg_welfare_weight_part_norm_C = stats.tot_welfare_weight_part_C/...
    (stats.tot_welfare_weight_part * stats.emp_C);

    % tax revenues
    stats.tax_rev_inc = derived_vars.tax_rev_inc;
    stats.exp_transfer = derived_vars.exp_transfer;
    stats.tax_rev_B = derived_vars.tax_rev_B;
    stats.tax_rev_E = derived_vars.tax_rev_E;
    stats.tax_rev_S = derived_vars.tax_rev_S;
    stats.tax_rev_total = derived_vars.tax_rev_total;
    stats.excess_revenue = derived_vars.excess_revenue;

    % percentiles
    handle_h = @(w) 1/stats.part .* econ.dist.h(w,stats.Y_M, ...
                                                stats.Y_R,stats.Y_C,par);
    pct = 0.1:0.1:0.9;
    [percentiles, nodes, weights] = ...
        tools.Integration.compute_percentiles(handle_h,pct,1,200,100);

    stats.part_M = econ.part_M(sol,par);
    stats.part_R = econ.part_R(sol,par);
    stats.part_C = econ.part_C(sol,par);

    % population shares
    integrand = @(w) econ.dist.dens_skill.f_M(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_M = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_R(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_R = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_C(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_C = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    percentiles_strct.percentile = pct';
    percentiles_strct.value = percentiles;
    percentiles_strct.q_B = econ.q_B * ones(size(pct'));
    percentiles_strct.q_E = econ.q_E * ones(size(pct'));
    percentiles_strct.q_S = econ.q_S * ones(size(pct'));
end